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The friction and adhesion between elastic bodies are strongly influenced by the roughness of the 
surfaces in contact. Here we develop a multiscale molecular dynamics approach to contact mechanics, 
which can be used also when the surfaces have roughness on many different length-scales, e.g., for 
self affine fractal surfaces. As an illustration we consider the contact between randomly rough 
surfaces, and show that the contact area varies linearly with the load for small load. We also 
analyze the contact morphology and the pressure distribution at different magnification, both with 
and without adhesion. The calculations are compared with analytical contact mechanics models 
based on continuum mechanics. 



1. Introduction 

Adhesion and friction between solid surfaces are com- 
mon phenomenons in nature and of extreme importance 
in biology and technology. Most surfaces of solids have 
roughness on many different length scales [3, Q, and it 
is usually necessary to consider many decades in length 
scale when describing the contact between solids 0. This 
makes it very hard to describe accurately the contact me- 
chanics between macroscopic solids using computer sim- 
ulation methods, e.g., atomistic molecular dynamics, or 
finite element calculations based on continuum mechan- 
ics. 

Consider a solid with a nominally flat surface. Let 
x, y, z be a coordinate system with the x, y plane parallel 
to the surface plane. Assume that z = fo(x) describe the 
surface height profile, where x = (x, y) is the position 
vector within the surface plane. The most important 
property characterizing a randomly rough surface is the 
surface roughness power spectrum C(q) defined byjT 



C(q) = -±pj#x (Mx)MO)}e* q x 
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Here (...) stands for ensemble average and we have as- 
sumed that h(s) is measured from the average surface 
plane so that (h) = 0. In what follows we will assume 
that the statistical properties of the surface are isotropic, 
in which case C(q) will only depend on the magnitude 
q = |q| of the wave vector q. 

Many surfaces tend to be nearly self-afflne fractal. A 
self-afflne fractal surface has the property that if part 
of the surface is magnified, with a magnification which 
in general is appropriately different in the perpendicu- 
lar direction to the surface as compared to the lateral 
directions, then the surface "looks the same", i.e., the 
statistical properties of the surface are invariant under 
the scale transformation^. For a self-afflne surface the 
power spectrum has the power-law behavior 
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FIG. 1: Surface roughness power spectrum of a surface which 
is self affine fractal for q\ > q > qo. The long-distance roll- 
off wave vector go and the short distance cut-off wave vector 
qi depend on the system under consideration. The slope of 
the logC — \ogq relation for q > qo determines the fractal 
exponent of the surface. The lateral size L of the surface (or 
of the studied surface region) determines the smallest possible 
wave vector q^ — 2tv/L. 



where the Hurst exponent H is related to the fractal di- 
mension D{ of the surface via H = 3 — D{. Of course, 
for real surfaces this relation only holds in some finite 
wave vector region qo < q < qi, and in a typical case 
C(q) has the form shown in Fig. ^ Note that in many 
cases there is a roll-off wavevector go below which C(q) 
is approximately constant. 

Let us consider the contact between an elastic solid 
with a flat surface and a hard randomly rough substrate. 
Fig. [21 shows the contact between the solids at increas- 
ing magnification £. At low magnification (£ = 1) it 
looks as if complete contact occurs between the solids at 
many macro asperity contact regions, but when the mag- 
nification is increased smaller length scale roughness is 
detected, and it is observed that only partial contact oc- 
curs at the asperities. In fact, if there would be no short 
distance cut-off the true contact area would vanish. In 
reality, however, a short distance cut-off will always exist 
since the shortest possible length is an atomic distance. 
In many cases the local pressure at asperity contact re- 
gions at high magnification will become so high that the 



material yields plastically before reaching the atomic di- 
mension. In these cases the size of the real contact area 
will be determined mainly by the yield stress of the solid. 

The magnification £ refers to some (arbitrary) chosen 
reference length scale. This could be, e.g., the lateral size 
L of the nominal contact area in which case £ = L/\, 
where A is the shortest wavelength roughness which can 
be resolved at magnification £. In this paper we will 
instead use the roll-off wavelength Ao = 27r/go as the 
reference length so that ( = Ao/A. 

Recently, a very general contact mechanics theory has 
been developed which can be applied to both stationary 
and sliding contact for viscoelastic solids (which includes 
elastic solids as a special case)jj. The theory was orig- 
inally developed in order to describe rubber friction on 
rough substrates. For elastic solids the theory can also 
be applied when the adhesional interaction is taken into 
account 0- In contrast to earlier contact mechanics the- 
ories, the theory presented in Ref. |4|, P is particularly 
accurate close to complete contact, as would be the case 
for, e.g., rubber on smooth surfaces. The basic idea be- 
hind the theory is to study the contact at different mag- 
nification. In particular, the theory describes the change 
in the stress distribution P(a, Q as the magnifications ( 
increases. Here 

P(a,C) = (S(a-a( X ,C))) (2) 

is the stress distribution at the interface when the sur- 
face roughness with wavelength smaller than A = Ao/C 
has been removed. In J2J), (...) stands for ensemble av- 
erage, and <j(x, () is the perpendicular stress at the in- 
terface when surface roughness with wavelength shorter 
than A = Ao/C nas been removed. It is clear that as the 
magnification ( increases, the distribution P(a, Q will be 
broader and broader and the theory describes this in de- 
tail. The (normalized) area of real contact (projected on 
the x?/-plane) at the magnification ( can be written as 

^1 = £daP(a,0- (3) 

where the lower integration limit + indicate that the 
delta function at the origin a = (arising from the non- 
contact area) should be excluded from the integral. The 
rubber friction theory described in Ref. depends on 
the function A(Q/Aq for all magnifications. This just re- 
flects the fact that the friction force results from the vis- 
coelastic deformations of the rubber on all length scales, 
and when evaluating the contribution to the friction from 
the viscoelastic deformations on the length scale A, it is 
necessary to know the contact between the rubber and 
the substrate at the magnification ( — Ao/A. Thus, not 
just the area of real (atomic) contact is of great inter- 
est, but many important applications require the whole 
function A((), and the pressure distribution P(cr, (). 

In order to accurately reproduce the contact mechan- 
ics between elastic solids, it is in general necessary to 
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FIG. 2: A rubber block (dotted area) in adhesive contact 
with a hard rough substrate (dashed area). The substrate 
has roughness on many different length scales and the rubber 
makes partial contact with the substrate on all length scales. 
When a contact area is studied at low magnification (£ = 1) 
it appears as if complete contact occurs in the macro asperity 
contact regions, but when the magnification is increased it is 
observed that in reality only partial contact occurs. 



consider solid blocks which extend a similar distance in 
the direction normal to the nominal contact area as the 
linear size of the contact area. This leads to an enor- 
mous number of atoms or dynamical variables already 
for relatively small systems. In this paper we develop 
a multiscale approach to contact mechanics where the 
number of dynamical variables scales like ~ N 2 rather 
than as ~ A 7 " 3 , where N x N is the number of atoms 
in the nominal contact area. As application we consider 
the contact mechanics between randomly rough surfaces 
both with and without adhesion, and compare the results 
with analytical contact mechanics theories. 

2. Multiscale molecular dynamics 

Let us discuss the minimum block-size necessary in a 
computer simulation for an accurate description of the 
contact mechanics between two semi-infinite elastic solids 
with nominal flat surfaces. Assume that the surface 
roughness power spectrum has a roll-off wavevector q = 
qo corresponding to the roll-off wavelength Ao = 27r/go- 
In this case the minimum block must extend L x « Ao 
and L y w Ao along the x and ^-directions. Furthermore, 
the block must extend at least a distance L z « Ao in 
the direction perpendicular to the nominal contact area. 
The latter follows from the fact that a periodic stress 
distribution with wavelength A acting on the surface of a 
semi-infinite elastic solid gives rise to a deformation field 
which extends a distance ~ A into the solid. Thus, the 
minimum block is a cube with the side L — Ao- 

As an example, if Ao corresponds to 1000 atomic spac- 
ings, one must at least consider a block with 1000 x 1000 
atoms within the xy-contact plane, i.e., one would need 
to study the elastic deformation in a cubic block with 
at least 10 9 atoms. However, it is possible to drasti- 
cally reduce the number of dynamical variables without 
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FIG. 3: Schematic structure of the model, (a) The fully 
atomistic model, (b) The multiscale smartblock model, where 
the solid in (a) is coarse grained by replacing groups of atoms 
with bigger "atoms". 

loss of accuracy if one notes that an interfacial roughness 
with wavelength A will give rise to a deformation field 
in the block which extends a distance A into the solid, 
and which varies spatially over a typical length scale A. 
Thus when we study the deformation a distance z into 
the block we do not need to describe the solid on the 
atomistic level, but we can coarse-grain the solid by re- 
placing groups of atoms with bigger "atoms" as indicated 
schematically in Fig. [3] If there are N x N atoms in the 
nominal contact area one need n ~ InTV "atomic" lay- 
ers in the z-direction. Moreover the number of atoms 
in each layer decreases in a geometric progression every 
time the coarse graining procedure is applied, so that the 
total number of particles is of order N 2 instead of N 3 . 
This results in a huge reduction in the computation time 
for large systems. This multiscale approach may be im- 
plemented in various ways, and in the Appendix A we 
outline the procedure we have used in this paper which 
we refer to as the smartblock. Another implementation 
similar to our approach can be found in Ref. Q. 

The smartblock model should accurately describe the 
deformations in the solids as long as the deformations 
varies slowly enough with time. However, the model can- 
not accurately describe the propagation of short wave- 
length phonons. This is, in fact, true with all forms 
of Hamiltonian multiscale descriptions of solids, because 
of the energy conservation together and the unavoid- 
able loss of information in the coarse grained region. In 
principle it should be possible to prevent the back re- 
flection of short wavelength phonons by describing the 
coarse grained region as a continuum, where the numer- 
ical calculation can be carried on through a Finite Ele- 
ment scheme. 0, S El [HI This indeed would require no 
coarse graining at all in the region treated with molecular 
dynamics, and a proper choice of the matching conditions 
between the atomistic and the continuum region. How- 
ever, with respect to contact mechanics and adhesion the 
back reflection of short wavelength phonons is not an im- 



portant limitation. With respect to sliding friction it may 
be a more severe limitation in some cases. 

Figure 03 illustrates a case where the block is in the 
form of a cube with atomically fiat surfaces. It is possi- 
ble to obtain curved surfaces of nearly arbitrary shape 
by "gluing" the upper surface of the block to a hard 
curved surface profile. This was described in detail in 
Ref. 4]. The elastic modulus and the shear modulus of 
the solid can be fixed at any value by proper choice of the 
elongation and bending spring constants for the springs 
between the atoms (see Ref. 4] and Appendix A). The 
upper surface of the smartblock can be moved with ar- 
bitrary velocity in any direction, or an external force of 
arbitrary magnitude can be applied to the upper surface 
of the smartblock. We have also studied sliding friction 
problems where the upper surface of the smartblock is 
connected to a spring which is pulled in some prescribed 
way. The computer code also allows for various lubricant 
fluids between the solid surfaces of the block and the 
substrate. Thus the present model is extremely flexible 
and can be used to study many interesting adhesion and 
friction phenomena, which we will report on elsewhere. 

We note that with respect to contact mechanics, when 
the slopes of the surfaces are small, i.e. when the surfaces 
are almost horizontal, one of the two surfaces can be con- 
sidered flat, while the profile of the other surface has to be 
replaced by the difference of the two original profiles [l]|. 
Thus, if the substrate has the profile z = fti(x) and the 
block has the profile z = /^(x), then we can replace 
the actual system with a fictive one where the block has 
an atomically smooth surface while the substrate profile 
h{x) = /i2(x) — fei(x). Furthermore, if the original solids 
have the elastic modulus E\ and E2 , and the Poisson ra- 
tio v\ and V2-, then the substrate in the fictive system 
can be treated as rigid and the block as elastic with the 
elastic modulus E and Poisson ratio v chosen so that 
(1 - p 2 )/E = (1 - vD/E, + (1 - vl)/E 2 . 

The results presented below have been obtained for 
a rigid and rough substrate. The atoms in the bottom 
layer of the block form a simple square lattice with lattice 
constant a. The lateral dimensions L x = N x a and L y = 
N y a. For the block, N x = 400 and N y = 400. Periodic 
boundary conditions are applied in the xy plane. The 
lateral size of the block is equal to that of substrate, but 
we use different lattice constant b « a/0, where <p = 
(1 + y/S)/2 is the golden mean, in order to avoid the 
formation of commensurate structures at the interface. 
The mass of a block atom is 197 a.m.u. and the lattice 
constant of the block is a = 2.6 A, reproducing the atomic 
mass and the density of gold. We consider solid blocks 
with two different Young's moduli: a hard solid with 
E = 77 GPa, like in gold, and a soft one with 0.5 GPa. 
The corresponding shear moduli were G = 27 GPa and 
0.18 GPa, respectively. 

The atoms at the interface between the block and the 




substrate interact with the potential 
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U(r) = 4e 
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where r is the distance between a pair of atoms. When 
a = 1, Eq. (|3J is the standard Lennard- Jones poten- 
The parameter e is the binding energy between two 
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atoms at the separation r = 2 1 / 6 ro- When we study con- 
tact mechanics without adhesion we put a = 0. In the 
calculations presented below we have used ro = 3.28 A 
and e = 18.6 meV, which (when a = 1 ) gives an inter- 
facial binding energy (per unit area) [l2j A7 « 4e/a 2 « 
11 meV/A 2 . 



3. Self affine fractal surfaces 

In our calculations we have used self affine fractal sur- 
faces generated as outlined in Ref. 3]. Thus, the surface 
height is written as 



ft(x) = ^£(q)e i [ q,x 



+0(q)] 



(5) 



where, since fe(x) is real, B(—q) = £>(q) and 0(— q) = 
— 0(q). If 0(q) are independent random variables, uni- 
formly distributed in the interval [0,27r[, then one can 
easily show that higher order correlation functions in- 
volving h(x) can be decomposed into a product of pair 
correlations, which implies that the height probability 
distribution = (S(h — h{x))) is Gaussian 3]. However, 
such surfaces can have arbitrary surface roughness power 
spectrum. To prove this, substitute (0 into QJ and use 
that 

( e i0(qO e i0(q")) = £ qW , 



gives 



C(q) = j^Jd 2 xJ2\B(<l')\ 
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where A is the nominal surface area, gives 

C(q) = ^l*(q)| 2 . 

Thus, if we choose 

B(q) = (2 7 r/L)[C(q)] 1 /2 ) 



(6) 



FIG. 4: (a) Fractal surface with the large wavevector cut-off 
qi = 2ir/b ~ 216 ^o- (b) The same surface as in (a) but at 
lower resolution with q± = 4qo. For a square 1040 A x 1040 A 
surface area. The fractal dimension Df — 2.2 and the root- 
mean-square roughness amplitude is 10 A. 



where L = A^ 2 , then the surface roughness profile (5J 
has the surface roughness power density C(q). If we as- 
sume that the statistical properties of the rough surface 
are isotropic, then C(q) = C(q) is a function of the mag- 
nitude q = |q|, but not of the direction of q. The ran- 
domly rough substrate surfaces used in our numerical cal- 
culations where generated using (0) and (0) and assuming 
that the surface roughness power spectra have the form 
shown in Fig. ^ with the fractal dimension D$ = 2.2 
and the roll-off wavevector qo — 3^l, where qL = 2tt/L x . 
We have chosen qo = 3qL rather than qo = qL since 
the former value gives some self- averaging and less noisy 
numerical results. We also used q\ = 2i\ jb w 216<2o (to- 
pography (a) in Fig. HJ) and some surfaces with several 
smaller values for q\ (Fig. 0] (b) shows the topography 
when qi = 4#o), corresponding to lower magnification 
(see Sec. 4). 

4. Numerical results 

In this section we illustrate our multiscale molecular 
dynamics (MD) approach by some applications. We first 
compare the MD results to two known contact mechan- 
ics results from continuum mechanics. Next we discuss 
contact mechanics for randomly rough surfaces both with 
and without adhesion. 

4.1. Test cases: Hertz contact and complete 
contact 

In 1881 Hertz presented an exact solution for the 
contact between two perfectly elastic solids with local 
quadratic profiles. The results were derived using the 
elastic continuum model and neglecting the adhesion be- 



tween the solids. In addition, Hertz assumed that the 
interfacial friction vanishes so that no shear stress can de- 
velop at the interface between the solids. When a spher- 
ical asperity is squeezed against a flat surface a circular 
contact area (radius t*h) is formed, where the pressure 
decreases continuously from the center r = to the pe- 
riphery r = rn of the contact according to 

r 2I 1 / 2 

Let us compare the prediction of our atomistic model 
with the Hertz theory. We use the Lennard- Jones poten- 
tial with a = 0, i.e. without the attractive term. In Fig. 
ISlwe compare the Hertz contact pressure (green line) with 
our numerical data (red data points). The numerical data 
were obtained for a rigid spherical tip squeezed against 
a flat elastic surface. Note that the pressure obtained 
from the MD calculation has a tail beyond the Hertz 
contact radius tr. Similar "pressure tails" were recently 
observed in molecular dynamics simulations by Luan and 
Robbins|]jj. The tail reflects the non-zero extent of the 
atom-atom interaction potential. The deviation between 
the molecular dynamics results and the continuum me- 
chanics results should decrease continuously as the size 
of the system increases. 

At the atomic level there is no unique way to define 
when two solids are in contact, and one may use sev- 
eral different criteria. One method is based on the force 
acting between the atoms at the interface and works best 
when the adhesional interaction is neglected. Thus, when 
two surfaces approach each other, the repulsive force be- 
tween the atoms increases. We may define contact when 
the repulsive force is above some critical value. When ad- 
hesion is included the interaction between the wall atoms 
becomes more long-ranged and it is less obvious how to 
define contact based on a force criterion, and we find it 
more convenient to use a criteria based on the nearest 
neighbor distance between atoms on the two surfaces. 
Thus, when the separation between two opposing surface 
atoms is less than some critical value, contact is defined 
to occur. However, we have found that neither of these 
two criteria gives fully satisfactory results. The reason 
is that if the critical force or the critical distance used 
to define when contact occurs is determined by fitting 
the Hertz pressure profile to the numerical data as in 
Fig. then the resulting values depend on the radius of 
curvature of the asperity. For example, for the Hertz con- 
tact in Fig. 0the contact area deduced from the atom- 
istic MD calculation agree with the Hertz theory if we 
choose the cut-off pressure p c « 0.7 GPa. However, if 
the radius of curvature of the asperity is 10 times smaller 
(R = 104 A) then, for the same penetration, the cut- 
off would be p c ~ 2.5 GPa, i.e., more than three times 
larger. On the other hand randomly rough surfaces have 
a wide distribution of curvatures and it is not clear how 
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FIG. 5: The pressure in the contact region between a spher- 
ical tip and a flat elastic surface. We show the simulation 
data and the theoretical Hertz result. The spherical tip has 
the radius of curvature R = 1040 A and the loading force 
4.6 x 10" 7 N. 

to choose the optimum cut-off distance or force. In this 
paper we have therefore used another way of determin- 
ing the contact area which turned out to be more unique. 
We will now describe this method. 

Let us consider the pressure distribution P(<r, Q at the 
interface. For Hertz contact we get the pressure distri- 
bution 

P(°) = T- I ( 8 ) 
A) Ja 

Using cr(x) from J7J) for r < th and cr(x) = for r > rn 
gives 

where A = i\r\ is the Hertz contact area. In Fig. El we 
show the pressure distribution in the contact region be- 
tween a hard spherical tip and an elastic solid with a flat 
surface. The red curve shows the simulation data, while 
the green curve is the theoretical Hertz result obtained 
by a suitable choice of A in Eq. Note that while 

the Hertz solution and the atomic MD simulation results 
agree very well for large pressure, there is a fundamental 
difference for small pressure. Thus, for the Hertz solu- 
tion, for small pressure a — > 0, P(<j) ~ a, while in the 
atomistic model P(c) increase monotonically as a — » 0. 
This difference is due to the long-range interaction be- 
tween the solid walls in the atomistic model, which is 
absent in the Hertz model. When the long range wall- 
wall interaction is taken into account the delta function 
at a = in the Hertz solution Q will broaden, resulting 
in a P(<j) which (for the small systems considered here) 
will decay monotonically with increasing a as observed 
for the atomistic model. Note that this effect is of exactly 
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FIG. 6: The pressure distribution in the contact region be- 
tween a spherical tip and a flat surface. We show the simula- 
tion data (red curves) and the theoretical Hertz result (green 
curves). Loading force in (a) is 4.6 x 10 -7 N and in (b) 
7.3 x 10" 7 N. 



the same origin as the "pressure tail" for r > rn in Fig. 

The fact that P(cr, () vanish linearly with a as a — > 
is an exact result in continuum mechanics with contact 
interaction (no long range wall- wall interaction), and is 
valid not just for the Hertz contact case, but holds in 
general [l^. However, as explained above, this effect will 
never be observed in the atomistic model if the wall-wall 
interaction is long-ranged. 

Note that the contact area A can be determined di- 
rectly by fitting the analytical expression for P(a) for 
the Hertz contact (Eq. (JJjJ)) to the numerical MD results 
for large enough pressures (see Fig. EJ- In the present 
case, for = 4.6 x 10 -7 N (Fig. EJa)) this gives a con- 
tact area A = i\r\ which is nearly identical to the one 
deduced from the fit in Fig. El A similar procedure can 
be used to determine the contact area between randomly 
rough surfaces using the following analytical expression 
derived from the contact mechanics theory of Persson 
(see Eq. (JT0|) below): 

p( a - 1 f p -(^o) 2 /4C7 _ p -(<r+ao) 2 /4CA 

1 ' g "2(rf)V2 r )> 

where cr is the nominal contact stress, and where the 
fitting parameter G = G(() can be related to the con- 
tact area using Eq. (3). Thus, if A/Aq « 1 we have 
G = (cjI/ii)(A/Aq)~ 2 . We have found (see below) that 
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FIG. 7: The normalized pressure distribution P(cr) at the 
interface between an elastic block (elastic modulus E = 
0.5 GPa) with a flat surface and a rigid randomly rough 
substrate. Because of adhesion complete contact occurs at 
the interface. The red curve is the simulation result and the 
green line is the Gaussian fit to the simulation data with the 
root-mean-square width cr rms = 0.229 GPa. The blue line is 
the theoretical Gaussian distribution obtained using contin- 
uum mechanics (see Appendix B). The theoretical rms width 
cr rms 0.164 GPa. 

this expression for P(a, Q can fit the numerical MD data 
very well (lending support for the accuracy of the Pers- 
son theory), and we have used this method to determine 
the contact area as a function of the squeezing force for 
randomly rough substrates. 

Let us consider the pressure distribution at the inter- 
face between a rigid randomly rough substrate and a fiat 
elastic surface when the solids are in complete contact. 
Complete contact can result either by squeezing the solids 
together by high enough force, or if the adhesional inter- 
action between the solids is high enough (or the elastic 
modulus small enough). However, when complete con- 
tact occurs the pressure distribution is the same. 

For an elastic solid with a flat surface in perfect contact 
with a hard randomly rough surface, continuum mechan- 
ics predict a Gaussian pressure distribution of the form 
(see Appendix B): 

p(rr) — - p -( cr - cr o) 2 /2cr r 2 ms 

n<7 > (2^)1/2^ e 

where the root-mean-square width <j rms is determined by 
the power spectrum: 

^rms = (<7 2 ) = 2 (1- V 2 ) 2 J ^ q3 °^ 

In Fig. E| we compare the theoretical pressure distribu- 
tion (blue curve) with the pressure distribution obtained 
from the atomistic model for the case where the com- 
plete contact results from the adhesive interaction be- 
tween the solids. The MD data are well fitted by a Gaus- 
sian curve, but the width of the curve is slightly larger 
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FIG. 8: Contact morphology for two different magnifications. 
The red color denotes contact regions for the low magnifica- 
tion C — 4, while the blue color corresponds to the contact 
regions for the high magnification £ = 216. 



than expected from the continuum mechanics theory 
cr rms (MD) = 0.229 GPa while a rms (theory) = 0.164 GPa. 
The randomly rough surface used in the MD calculation 
is self affine fractal the whole way down to the atomic dis- 
tance, and one can therefore not expect the continuum 
mechanics result for P(cr), which assumes "smooth" sur- 
face roughness, to agree perfectly with the MD result. 

4.2. Contact mechanics without adhesion 

Here we study contact mechanics without adhesion as 
obtained with a = in Eq. (4), corresponding to purely 
repulsive interaction between the walls. Fig. [HI shows the 
contact morphologies at different magnifications £ for the 
same load. The red and blue color indicate the contact 
area at low (£ = 4) and high (£ = 216) magnification, re- 
spectively. Note that with increasing magnification the 
contact area decreases, and the boundary line of the con- 
tact islands becomes roue; her. In Ref. Q and it has 
been shown that the statistical properties of the contact 
regions exhibit power-law scaling behavior. At low mag- 
nification ((" = 4) it looks as if complete contact occurs 
between the solids at asperity contact regions. However, 
when the magnification is increased, smaller length scale 
roughness is detected and it is observed that only partial 
contact occurs at the asperities. In fact, if there were no 
short distance cut-off in the surface roughness, the true 
contact area would eventually vanish. But in reality a 
short distance cut-off always exists, e.g. the interatomic 
distance. 

Fig. [HI shows the pressure distribution in the contact 
area for two different magnifications. When we study 
contact on shorter and shorter length scale, which corre- 
sponds to increasing magnification the pressure distri- 
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FIG. 9: The pressure distribution in the contact area for 
two different magnifications. The red line corresponds to the 
pressure distribution for low magnification £ = 4, while the 
green line is for high magnification £ = 216. 
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FIG. 10: The relative contact area A/Ao, as a function of 
applied stress Fn/^o- Results are presented for two different 
magnifications £ = Ao/A = 4 and 32. The fractal dimension 
isD/ = 2.2. 



bution becomes broader and broader. 

Fig. El shows that the contact area varies (approxi- 
mately) linearly with the load for the small load at two 
different magnifications £ = 4 and 32. The contact area 
was determined as described in Sec. 4.1. by fitting the 
pressure distribution to a function of the form ([TUj). The 
pressure distributions and the fitting functions are shown 
in Fig. ITT1 and fT21 for £ = 4 and 32, respectively. The 
slope of the lines in Fig. Ellis on ly a factor 1.14 larger 
than predicted by the contact theory of Persson (see Sec. 
5). 

In Fig. ^|we show the variation of the contact area 
with the nominal squeezing pressure for the highest mag- 
nification case £ = 216. In this case we have defined con- 
tact to occur when the separation between the surfaces 
is below some critical value r c = 4.3615 A. In contrast to 
the definition used above, this definition does not give a 
strict linear dependence of the contact area on the load 
for small load as found above when the contact area is 
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FIG. 11: The stress distribution for f = 4 for three different 
nominal pressure. 
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FIG. 12: The stress distribution for £ = 32 for three different 
nominal pressure. 
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FIG. 13: The relative contact area A/Ao, as a function of 
applied stress F^/Ao. Results are presented for the highest 
magnification £ = 216. Contact is defined when the separa- 
tion between the surfaces is below a critical value. The fractal 
dimension is Df — 2.2. 
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FIG. 14: Contact morphology with adhesion and without 
adhesion. The blue color region denotes the contact without 
adhesion. The red color region denote the additional contact 
area when the adhesional interaction is included. 



defined using the stress distribution. 

4.3. Contact mechanics with adhesion 

In this section we include the adhesive interaction i.e. 
we put a = 1 in Eq. (4). Fig. ^] presents the contact 
morphology both with and without the adhesion at the 
highest magnification (£ = 216). The regions with blue 
color denotes the contact area without adhesion. The 
red color region denotes the additional contact area when 
adhesion is included. The contact area with adhesion is, 
of course, larger than that without adhesion since the 
attractive adhesional interaction will effectively increase 
the external load[I3, [H, 13 - 

Fig. [131 shows the pressure distribution P(cr, () at high 
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FIG. 15: The pressure distribution with and without adhe- 
sion. The red curve denotes the pressure distribution with 
adhesion while the green curve is without adhesion. 



magnification with and without adhesion. When adhe- 
sion is neglected (corresponding to the a = in (4)), the 
pressure is positive in the contact area and P(a, £) = 
for a < 0. When the adhesive interaction is included, the 
stress becomes tensile close to the edges of every contact 
region and P(cr, () is in general finite also for a < 0. 

5. Discussion 

Several analytical theories, based on continuum me- 
chanics, have been developed to describe the contact 
between elastic bodies both with and without the ad- 
hesional interaction. Here we will compare the results 
presented above with the predictions of some of these 
theories. 

Persson has developed a contact mechanics theory 
where the surfaces are studied at different magnification 
C = ^o/A, where Ao is the roll-off wavelength and A the 
shortest wavelength roughness which can be observed at 
the magnification £. In this theory 4] the stress distri- 
bution P(a, C) at the interface between the block and 
the substrate has been shown to obey (approximately) 
a diffusion-like equation where time is replaced by mag- 
nification and spatial coordinate by the stress a. When 
the magnification is so small that no atomic structure 
can be detected, the surface roughness will be smooth 
(no abrupt or step-like changes in the height profile) and 
one can then showfl^ that in the absence of adhesion 
P(0, C) = 0. Using this boundary condition the solution 
to the diffusion-like equation gives the pressure distribu- 
tion at the interface (a > 0): 
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Substituting into (|12|) gives after some simplifica- 
tions 
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Thus, for small nominal squeezing pressure cfq <C G 1//2 
we get 



(14) 



Ao (ttG) 1 / 2 ' 
Since the squeezing force = (JqAq we can also write 

-1/2 



(15) 



(10) 



where E* = E/{1 - v 2 ) and k = (8/tt) 1 / 2 . Thus, for 
small squeezing force Fn the theory predicts a linear de- 
pendence of the area of real contact on the load. 

For very high squeezing force <Jo ^> G 1 ^ 2 complete con- 
tact will occur at the interface. In this case the second 
term on the rhs in (fTU|) can be neglected, so the pressure 
distribution is a Gaussian centered at ao and with the 
root-mean-square width <r rms = (2G) 1/2 . This result is 
exact (see Appendix B). Thus, the theory of Persson is 
expected to give a good description of the contact me- 
chanics for all squeezing forces. All other analytical con- 
tact mechanics theories are only valid when the squeezing 
force is so small that the area of real contact is (nearly) 
proportional to F^. But in many important applications, 
e.g., in the context of rubber friction and rubber adhe- 
sion, the area of real contact for smooth surfaces is often 
close to the nominal contact area. 

The standard theory of Greenwood and Williamson 
[2fl| describe the contact between rough surfaces (in the 
absence of adhesion), where the asperities are approxi- 
mated by spherical cups with equal radius of curvature 
but with Gaussian distributed heights. In this theory 
the area of real contact dependent (slightly) non-linearly 
on the load for small load, and can therefore not be di- 
rectly compared with the Persson result ([T5|). Bush et al 
[2~H developed a more general and accurate contact the- 
ory. They assumed that the rough surface consists of a 
mean plane with hills and valleys randomly distributed 
on it. The summits of these hills are approximated by 
paraboloids, whose distributions and principal curvatures 
are obtained from the random precess theory. As a result 



of more random nature of the surface, Bush et al found 
that at small load the area of contact depends linearly on 
the load according to (|T5|) but with k = (2tt) 1 ^ 2 . Thus 
the contact area of Persson's theory is a factor of 2/ir 
smaller than that predicted by Bush. Both the theory 
of Greenwood and Williamson and the theory of Bush 
et al assume that the asperity contact regions are inde- 
pendent. However, as discussed in [l4|, for real surfaces 
(which always have surface roughness on many different 
length scales) this will never be the case even at a very 
low nominal contact pressure, which may be the origin 
of difference of 2/ir between Persson's theory and Bush's 
theory. 

Hyun et al performed a finite-element analysis of con- 
tact between elastic self-afflne fractal surfaces 16]. The 
simulations were done for rough elastic surface contact- 
ing a perfectly rigid flat surface. They found that the 
contact area varies linearly with the load for small load. 
The factor n was found to be between the results of the 
Bush and Persson theories for all fractal dimensions Df . 
For Df = 2.2 (corresponding to H = 0.8) they found that 
n was only ~ 13% larger than predicted by the Persson 
theory. 

The red curves in Fig. EH shows the pressure distri- 
bution from the simulations for several different values 
of the magnification £ = qi/q$ = 4, 8, 32 and 216, ne- 
glecting the adhesion. In the simulations the nominal 
squeezing pressure ao = 800 MPa. The best fit (green 
curves in Fig. IT6|) of the pressure distribution (fTUj) to the 
numerical results is obtained if G -1 / 2 is taken to be a 
factor 1.14 larger than predicted by the Persson theory 
[Eq. (10)], corresponding to a contact area which is 14% 
larger than predicted by the analytical theory, in good 
agreement with the results obtained by Hyun et al. 

Our simulations show that the contact area varies lin- 
early with the load for small load, see Fig. Figs. ITUl 
andlT?)lshow that the slope a(() of the line A = a(()F de- 
creases with increasing magnification as predicted by 
the analytical theory QEl- Thus, while A/A = 0.072 
for C = 4 we get A/A = 0.038 for ( = 32, which both 
are 14% larger than predicted by Eq. (p"3|) . 

6. Summary and conclusion 

In this paper we have developed a Molecular Dynamics 
multiscale model, which we have used to study the con- 
tact between surfaces which are rough on many different 
length scales. We have studied the contact morphologies 
both at high and low magnification, with and without 
adhesion. We have shown that in atomistic models it is 
a non-trivial problem how to define the area of real con- 
tact between two solids. Our study shows that the area 
of real contact is best defined by studying the interfacial 
pressure distribution, and fitting it to an analytical ex- 
pression. The numerical results are consistent with the 
theoretical results that the contact area varies linearly 




^=4 

A/A o=0.072 



Q_ 
Q_ 



-13 



C=8 

A/A o=0.057 




IHU. I I 



^=32 

A/A o=0.038 



13 



-9 



o- -11 



-13 



C=216 

A/A o=0.025 



20 40 60 80 

pressure (GPa) 

FIG. 16: The pressure distribution at four different magnifi- 
cations C — Qi/Qo = 4, 8, 32 and 216 for the squeezing pressure 
ao = 800 MPa. The red curves is the pressure distribution ob- 
tained from the computer simulation, while the green curves 
is from the analytical theory assuming that G~ 1 ^ 2 , and hence 
the relative contact area, is a factor of 1.14 larger than pre- 
dicted by the analytical theory, Eq. (10). 



with the load for small load, where the proportionality 
constant depends on the magnification L/X. For a ran- 
domly rough surfaces with the fractal dimension Df = 2.2 
(which is typical for many real surfaces, e.g., produced 
by fracture or by blasting with small particles) we have 
found that for small load (where the contact area is pro- 
portional to the load) the numerical study gives an area 
of atomic contact which is only ~ 14% larger than pre- 
dicted by the analytical theory of Persson. Since the 
Persson's theory is exact in the limit of complete con- 
tact, it is likely that the Persson theory is even better for 
higher squeezing loads. 
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FIG. 17: The model of Persson and Ballone with long range 
elasticity. Side view. 
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Appendix A: The "smartblock" model for mul- 
tiscale molecular dynamics 

Here we present a detailed description of the multiscale 
model implemented in our Molecular Dynamics (MD) 
simulations. Persson and Ballone 22] introduced a sim- 
ple and effective model to study the boundary lubrication 
between elastic walls. For each wall only the outermost 
layer of atoms was considered. These atoms were able to 
interact with the lubricant or with the atoms of the other 
wall with Lennard- Jones potentials. The walls' atoms 
were connected to a rigid surface through special springs 
which exert an elastic reaction not only to elongation, 
but also to lateral bending. The walls' atoms are cou- 
pled with their in-plane neighbors with similar springs. 
It is also possible to use curved elastic walls by connect- 
ing the vertical springs to a curved rigid surface rather 
than a flat surface as in Fig. El 

The model of Persson and Ballone catches two essen- 
tial features: Firstly the walls are not rigid, they can de- 
form (differently from previous models) and the descrip- 
tion takes into account the elastic energy stored during 
compression or stretching, which is an essential ingredi- 
ent for the study of the squeeze-out. Secondly both the 
Young modulus and the shear modulus can be indepen- 
dently tuned via the choice of the elastic constants of the 
springs. 

The model of Persson and Ballone works well when the 
solid is exposed to uniform shear or uniform elongation 
or compression. However, when there are spatial varia- 
tions in the stress at the surface of the solid, for instance 
when the displacement at the interface comprises short 
wavelength Fourier components, then the model does not 
allow a proper description of the elastic deformation field. 
In particular, when a periodic stress acts on the surface 
of an elastic solid, the displacement field decays expo- 



nentially into the solid, and this aspect is absent in the 
Persson-Ballone model. 

The solution to overcome this limitation is straight- 
forward: we explicitly introduce many layers of atoms, 
placed on the points of a simple cubic lattice, and cou- 
pled with springs to their nearest neighbors. 

The "springs", as in the previous model, are special, 
since they can resist to lateral bending. The force due to 
a vertical spring connecting two consecutive atoms 1 and 
2 along the z axis is given by the formulas below, where 
a is the lattice spacing, that is the equilibrium length of 
the spring: 

F x = -k b Ax = -k b (x 2 - xi), (Al) 
F y = -k b Ay = -k b (y 2 - j/i), (A2) 

F z = -k(Az -a) = -k[z 2 - (z ± + a)]. (A3) 

Analogous formulas hold for the springs parallel to the y 
and to the z axes. The two elastic constants of the spring, 
namely k and k b: are related to the Young modulus E and 
the shear modulus G respectively: k = Ea and k b = Ga. 

In some circumstances it is useful to simulate quite 
large and thick samples. Moreover high resolution up to 
the atomic level is needed in part of the sample, typically 
at the interface. The solution to avoid excessive compu- 
tational time is a multiscale approach: high resolution is 
achieved where it is needed, but a coarse grained descrip- 
tion is employed when it is feasible. The coarse graining 
can happen more times, and to various degrees of res- 
olution, so that a multilevel description of the system 
comprising many hierarchies is implemented. 

The grid structure of the smartblock allows a simple 
procedure to achieve a multiscale description: groups of 
atoms can be replaced by single, bigger atoms, and the 
elastic constants of the springs are redefined to guaran- 
tee the same elastic response. In many calculations per- 
formed by our group we used to replace a cube of 2 x 2 x 2 
particles with a single particle, repeating this merging 
procedure every two layers. More generally any change 
of resolution involves merging together a box made of 
m x x m y x m z particles. The three numbers m X: m y and 
m z are called merging factors along the three axes. 

The equilibrium position of the new particle is in the 
center of mass of the group of particles merged together. 
Its mass is m x m y m z times the mass of the original par- 
ticles, so that the density does not change. In fact the 
masses are only important to study the kinetic, but they 
do not influence the static equilibrium configuration. 

The three merging factors can be chosen indepen- 
dently. The easiest way to calculate the new springs' 
elastic constants is by considering the merging only along 
one of the axes. Fig. ^| sketches the case m z = 2, 
m x = m y = 1 (no change of lattice constant along x 



FIG. 19: Change of lattice spacing along a direction parallel 
to the interface between the two grids. 



FIG. 18: The grid of particles is coarse grained by replacing 
two atoms with a single one. Merging factors m x = m y — 1, 
m z = 2. Masses, equilibrium positions and spring constants 
are changed accordingly. 



and y). Along the direction of merging the new spring 
constants for elongation and bending are k' = k/m z and 
k' h = kb/m z respectively. The longer springs get pro- 
portionally smaller elastic constants, as it happens when 
springs are connected in series. In the two directions or- 
thogonal each spring replaces m z old springs in parallel 
configuration, so the elastic constants increase propor- 
tionally: k' = m x k, k' h = m x kb. Below there is the 
general formula giving the new elastic constants of the 
springs along the z axis, with arbitrary merging factors: 



k' 



m x m y . 



k' h 



m x m y 



(A4) 



Analogous formulas hold for the springs parallel to the x 
and y axes. 

To get the whole picture we have to characterize the 
springs at the interface between the two lattices, e.g., 
the ones crossing the dashed line in Fig. ^1 When the 
merging is in the direction z orthogonal to the interface 
both elastic constants k and k^ get multiplied by the 
factor 2/(1 -\-m z ). Actually their length is \{m z + l)a z , 
a z being the old lattice constant along z. Each of these 
interface springs can be thought as half a spring of the 
old grid connected with half a spring of the new grid. 

When the merging is along a direction orthogonal 
to the interface between the two grids, as sketched in 
Fig. [Hll then the spring constants do not change, but 



the forces between the particles are calculated taking 
into account the in-plane shift between the atoms of 
the two grids. Each interface particle of the upper lat- 
tice interacts with m x x m y particles of the lower lat- 
tice. The equations (Al) and (A2) are modified: F x — 
—kb(Ax + x-shift), F y = — kb(Ax + y-shift). The two in- 
plane shifts depend on the pair of particles considered. 

Appendix B: Pressure distribution at complete 
contact between randomly rough surfaces 

Here we calculate the pressure distribution at the inter- 
face between two solids in complete contact. We assume 
that one solid is rigid and randomly rough and the other 
solid elastic with a flat surface. The pressure distribution 
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where 



F(a) = (e 



where cr(x) is the fluctuating pressure at the interface. 
Next, writing 



cr(x) = J d 2 q <r(q)e* q 
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where we have used the relation between <r(q) and the 
Fourier transform h(q) of the height profile h(x) derived 
in Ref. |4(, we get 

F = ^exp (-ia J d 2 q -J^h^e** 



Next, using that h(q) are independent random variables 
we get 



F = e -« 2 e/2 



where 
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However (see Ref. j^) 

(h(q)h(q'))=C(q)5(q + q') 

so that 
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Substituting (B2) in (Bl) and performing the a-integral 
and using (B3) gives 



P(a) 
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where the root-mean-square width <j rms is determined by 
the power spectrum: 

°L* = = 2 (1-1/2)2 J <? C ^) 
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